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Abstract. An optimal algorithm is described for solving the deconvolution 
problem of the form kw :— k{t — s)u{s)ds = f(t) given the noisy data fs, \\f — 
fs\ \ ^ S. The idea of the method consists of the representation k = A{I + S), where 
S' is a compact operator, / -t- S* is injective, / is the identity operator, A is not 
boundedly invertible, and an optimal regularizer is constructed for A. The optimal 
regularizer is constructed using the results of the paper MR 40#5130. 

1. Introduction 

Deconvolution problem consists of solving equation of the form 

(1.1) ku:= [ k{t-s)u{s)ds:=ki.u^ f{t), 0<t<T, 

Jo 

where k{t), i > 0, is a kernel of linear integral equation Hl.l|l . k-ku is the convolution. 
It is important in many engineering applications, in physics, and other areas. There 
is a vast literature on deconvolution methods, see, for example, 6 . 

If the operator k in ^1.1^ is considered as an operator on X :— L°°{0,T), and 
/p \k{t) \ dt < oo, then k is not boundedly invertible, so problem l|l.HI is ill-posed. 
Assume that the data / are noisy: fs is given, such that \\f — fs\\ < S. In this 
case it is natural to seek an approximate solution of equation (|l.ll) in the class 
Qs {u E X : \\ku — fsW < S}. However, for ill-posed equation (|1.1() an arbitrary 
element us € Qs cannot be taken as an approximate solution to (|l.l|l . since us is 
not continuous with respect to S in general. In order to select possible solutions 
one needs to use a priori information (usually available) about the solution, which 
may be of a quantitative or qualitative nature. 

The usage of qualitative a priori information makes it possible to narrow the 
class of solutions, for example, to a compact set, so that the problem becomes stable 
under small changes in the data. This leads to a concept of a quasisolution |H| . 
Various algorithms for approximate determination of quasisolutions were studied 
inlHl. 

1 



2 



A priori information of a qualitative nature (for example, smoothness of the 
solution) generates different approaches. The one which is used often is varia- 
tional regularization [20], JO], which allows one to construct stable approximate 
solutions to ill-posed problems by means of a stabilizing functional. The varia- 
tional method has been extensively developed in and certain a priori and 
a posteriori choices of a regularization parameter e = e{S) have been designed and 
implemented |ni,|2I- 

One can also find approximate solutions to by iterations (see [2], PP), 

taking Xn — R{fs, Xn-i, Xn-k), where k < n. For these solutions to be stable 
under small changes of the data, the iteration number n — n{d) yielding x„ must 
depend on the 6 suitably. 

Other important techniques in theory of ill-posed problems give regularizing op- 
erators by using Fourier, Laplace, Mellin, and other integral transforms, statistical 
regularization, and the dynamical systems method (DSM) 

In ri4 some general new approaches are proposed for solving an ill-posed decon- 
volution problem. One of these approaches is based on the following idea. Assume 
that the operator k in 1)1. l|l can be decomposed into a sum k := A -t- B, where the 
operator A~^B := S is compact in the Banach space X, in which k acts, and I S 
is boundedly invertible. By the Fredholm alternative, it is equivalent to assuming 
that N{I + S) — {0}, where M{A) is the null space of A. In this case / + S* is an 
isomorphism of X onto X, 'R.{A) = TZ(k), where 7^(A) is the range of the operator 
A, and 

(1.2) ku = A{I + S)u^ fs. 

If a regularizer for A is known, then H1.2|l can be solved stably by the scheme 

(1.3) us^{I + S)-^R{6)fs, 
and 

(1.4) \\u-us\\^0 as S^O. 

Since / + 5 is an isomorphism, the error ||t; — uajl of the approximation of the 
solution of the equation Av — fs by the formula vs — R{S)fs is of the same order 
as — In this paper (see sections 2 and 3) we show that the proposed method 
is practically efficient and works better than the variational regularization. 

Theoretically the proposed method is optimal on the class of the data defined as 
a triple {S,fs,M2}, where / G C^(0,T), ||/"|| < M2, and / is otherwise arbitrary, 
fs £ i°°(0,T) and ||/ — /^H < S and fs is otherwise arbitrary. 

The operator R{6), defined in (2.3) and originally proposed in [10] for stable 
numerical differentiation, yields an optimal estimate of /' in L°°(0, T)— norm in 
the following sense: 

inf sup \\Tfs- f'\\>{2M2S)^/^ 

^ {fs--\\f-fs\\<S. ||/||<A/2} 

where the infimum is taken over all, linear and non-linear, operators T : X — > X, 
X = L°°(0,T), the supremum is taken over all / and fs subject to the conditions 
f£C\0,T), liril < Af2, \\f~fs\\<S, and 

\\R{S)fs-f'\\<{2M2S)'/\ 

(see e.g., [16],[17], [14]). 
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This argument shows that our " deconvohition" method for stable solution of 
(1.1) is optimal on the above data set: the operator R{S) gives an optimal (on the 
above data set) approximation of /'. Inversion of an isomorphism I + S, where S 
is a compact operator, can be done very accurately by a projection method, for 
example, so that the total error of the solution is of the same order as the error 
obtained by applying R{6). 

2. The case k{t) e C^{0,T) 



Let k{t) G C^(0,T) and k{0) ^ 0. Then without loss of generality one can take 
k{0) ==1. As in .14,, write (|l.l(l as 

(2.1) ku= [ u{s)ds+ [ [k{t- s) -l]u{s) -.^ Au + Bu^ f. 

Jo Jq 

Assume that f{x) is given by its (5-approximation, i.e. one knows fsix) such that 
1/ — fs\\x < S. In the experiments of this section S = 0.1. Let A~^B := S. Then 

(2.2) ku^ A{I + S)u^ f. 

Stable inversion of A is equivalent to stable numerical differentiation of noisy data, 
and therefore as a regularizer R{S)fs for A one can use (see 53 1 El i CHI: QHIi and 
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also 

fsit + hiS))-fs{t-hi6)) 

(2.3) R{S)fs . 

with h{6) = {^y^\ < ^2. Hence 

(2.4) iI + S)us^R{S)fs, 



where S* is a Volterra operator: Sug — J^k'{t — s)us{s) ds. To test numerical 
efficiency of the above deconvolution algorithm, we take 

(2.5) fc(j/) = exp(ay), fjt) = + "^^^'P^^^^ " ^^'^ + ' ""^^^^ . 

Then equation has the exact solution: 

(2.6) Uorig{t) = sin(6f) + cos(&t). 

The graphs of / and its (5-approximation, /i, for T — \, a — \, h ~ 27r, are pre- 
sented in Figure 1. The perturbation was generated as a sum of five sinusoids with 
various periods and amplitudes in such a way that ||/ — fs\\x < 0.1. For 5 = 0.1 

/ \ 1/2 

and for the above choice of /, T, a, and b, one has h{S) ^ [jj^j — 0.1253. Since 
in practice often only an estimate for M2 may be available, our first experiment was 
done with the approximate value of h(S), namely h = 0.105. The goal of the first 
experiment was to compare the results obtained by the deconvolution method sug- 
gested in jl4 | and by the variational regularization with a choice of the parameter 
by the Morozov discrepancy principle. The integral in Hl.l|l was calculated by the 
corrected trapezoid formula (see [7]) with the number of node points n — 200 on 
the interval [0, 1]. The graphs of Ujjis(;(i) and ^^^^^-^^{t) as well as the graph of the 
original solution, MQj.jg(t), for h{d) = 0.105 and n = 200 are given in Figure 2. One 
can see from the picture that method 14 provides higher quality of reconstruction. 



Table 1. 



t 


""exact (^) 


"disc(*) 


'"deconv(*) 


0.05 


1.26007351067010 


0.88613219081253 


1.61104047434242 


0.15 


1.39680224666742 


0.77345683250358 


1.16771020714854 


0.25 


1.00000000000000 


0.78531546607804 


0.97567292365993 


0.35 


0.22123174208247 


0.32264819143761 


0.46901890046136 


0.45 


-0.64203952192021 


-0.01522580641369 


-0.94010917284100 


0.55 


-1.26007351067010 


-0.72058597578420 


-1.39931254313538 


0.65 


-1.39680224666742 


-0.65363525334725 


-1.13246945274454 


0.75 


-1.00000000000000 


-0.84181827797783 


-1.26012127085008 


0.85 


-0.22123174208247 


-0.48659287989254 


-0.24854842261471 


0.95 


0.64203952192021 


-0.25478764331776 


0.99713489435843 
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Figure 2 
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h=0.105, n=200 

Table 1 allows one to analyze the computed values of '"(ijs(;;(i) and U(jg(,Q-^y(t) for 
h{S) = 0.1 and n = 10. The regularizatfon parameter for the variatfonal regulariza- 
tion calculated by the Morozov discrepancy principle, e^isc equal to 0.0275 for 
our particular fg. The functions U(\[g(.{t) and ""dcconv^'') approximate the exact 
solution UQj-ig(i) with the relative errors S^[g(, — 0.5216 and <5(jg(jonv ~ 0.2470, 
respectively, for n — 10. 



Table 2. 



n 


'^disc 


deconv 


10 


0.52160739359373 


0.24703402714545 


50 


0.47882066139400 


0.29886579484582 


100 


0.49421933901812 


0.29887532874922 



The deconvolution procedure j21 is applicable when the constant Ma, a > 1 is 
known. Here Ma is the bound on the /'•°\ a > is a real number, and /'•°-' is 
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Figure 6 





the (fractional order) derivative of / (see for details). Figures 3-6 show the 
dependence of the quality of calculations provided by the deconvolution technique 
for different values of h{d) with the same fs that is given in Figure 1. The level of 
reconstruction is acceptable for all values of h{S) £ (0.09, 0.3), but the best quahty 
is attained for the near-optimal values: h{S) = 0.1 and h{S) = 0.2. Outside the 
interval (0.09, 0.3) the reconstruction by the variational regularization works better 
because h{S) is away from its optimal value. 

Table 2 contains relative errors, ^^^jgc a-nd ^^ieconv' values of n = 10, 50, 100. 
In both cases the relative errors are not decaying further as n increases, because 
the major component in these errors come from the noise level, and not from the 
error of the computational methods. 



3. Kernel of the type k{t) = + m{t), < 7 < 1, m{t) e 



In this section we solve (|l.l(l with the kernel k{t) of the form 



(3.1) 



r(7) 



As in ,14,, write Hl.l|l as 
(3.2) 



+ m(t), 0<7<1, m(t) e 



kw :— Au + Bu ~ /, 
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Figure 7 
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Figure 8 




0.4 0.5 
delta=0.1 



where 
(3.3) 



^7-1 

Au :— „ , , * It, Bu := m-k u. 



One has ([S], pp. 117-118) A ~ r{i-j) So (t-s)-' ^™ce the right-hand side / 
is given by its (5-approximation fs, \\f — fs\\x < S,we replace A~ ^ by the regularizer 
RiiS) (see d): 



(3.4) 



Ri{S)fs :-- 



1 



■ ds. 



r(i - 7) J, it - s)-' 

The operator R{S) in (|?^ is defined by formula with /i = 0.12. One gets 

(3.5) iI + S)us = Ri{S)fs. 

and 

.ofiN c 4-1 D 1 m{0)us{s)+ J^m'{s~p)us{p)dp 
3.6 5w5 A Bus = / jf ds. 

r(i-7)yo [t-sp 

The goal of the experiment was to compare two numerical methods for solving 
(ll.ll) - (|3.1|) : deconvolution method H3.5|) - (|3.t)|) and variational regularization with a 
choice of the parameter by the discrepancy principle. 
The function 



^^^> " WTT) V (l + 7)(2 + 7); ^ 3" " loj 



i G [0,1], 
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was chosen as the solution to direct problem Hl.l|l - H3.1|l with m{t) = and the 
model function Uexact{t) = 1 — i^. Then for the numerical tests the noisy function 
fs, \\f — fs\\x < S, 6 = 0.1, was used. The graphs of / and fs for 7 = 0.1 are given 
in Figures 7 and 8. 



Table 3. 



t 


""exact (^) 


"disc(*) 


'"deconv(*) 


0.05 


0.99000000000000 


0.93361656127658 


1.00281244943820 


0.15 


0.96000000000000 


0.85983757148008 


0.98540317766041 


0.25 


0.91000000000000 


0.78772680932067 


0.93442776030698 


0.35 


0.84000000000000 


0.71362820985483 


0.85896131974899 


0.45 


0.75000000000000 


0.63504318796224 


0.76170936024357 


0.55 


0.64000000000000 


0.55028612377340 


0.64396777696250 


0.65 


0.51000000000000 


0.45824705747747 


0.50654835394340 


0.75 


0.36000000000000 


0.35823345537370 


0.35004574315540 


0.85 


0.19000000000000 


0.24979168725846 


0.17493299888351 


0.95 





0.13214536398250 


-0.01840021170081 
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Figures 10-12 illustrate the numerical performance of method H3.5() - (|3.6() and 
variational regularization for 7 = 0.1. The solutions evaluated by formulas (|3.5() - 
H3.6(l and by variational regularization for n = 10 and 7 = 0.1 are also presented in 
Table 3. 

The results obtained for our particular test problem show that for small values 
of 7 the deconvolution approach is superior to variational regularization both in 
terms of accuracy and stability. However as 7 is getting bigger, the efficiency of 
the deconvolution method (as well as the efficiency of variational regularization) is 
getting worse. This is happening because when 7 is close to 1, the ill-posedness of 
problem Hl.l|l - H3.1|l grows due to the errors in calculations of the singular integral. 
One can compare Figures 9-12 and 13-16. Moreover, as 7 changes from 0.1 to 0.9, 
method (|3.5|I - H3.6|I becomes very sensitive to slight variations of h{S). To illustrate 
this phenomena, we present the dependence of relative errors and discrepancies on 
h{d) for 7 = 0.1 and 7 = 0.5 in Figures 17 and 18. For 7 = 0.1 the relative error 
of the deconvolution method remains less than 10% when h{5) e (0.05,0.12), while 
for g = 0.5 the relative error is only small for h = 0.1. 

Finally, it is important to mention that CPU time for both methods, H3.5|) - (|3.t)|) 
and variational regularization, is approximately the same and it is very small: about 
3-4 milliseconds for n = 200. 

Conclusion. The paper presents numerical results of the implementation of the 
deconvolution method developed by AGR and presented together with other results 
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Figure 1 7 




0.05 0.1 0.15 0.2 0.25 

Deconvolution method for gamma=0.5: errors vs. h(delta) 

in [14]. The method is shown to be optimal in the sense explained in Section 1. 
The numerical results confirm the theoretical results on which the method is based. 
It is shown that the method is more accurate than the variational regularization 
method with the regularization parameter chosen by the discrepancy method. 
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